En este ejercicio vamos a estudiar la tasa de enfermedades transmitidas por los alimentos causadas por el consumo de hortalizas y frutas en la Unión Europea y Estados Unidos con los datos del trabajo de Callejón et al.[1] (2015).
Las enfermedades del análisis fueron Norovirus, Salmonella spp, Escherichia coli, Campylobacter spp, Shigella spp, Clostridium spp, Staphylococcus spp, Yersinia spp, Bacillus spp, otros virus y otros microorganismos. Las verduras eran Ensalada (todos los productos relacionados con la ensalada), Hojas (todos los productos relacionados con las hojas), Tomate y Otras verduras. Las frutas fueron Germinados (todos los productos relacionados con los germinados), Bayas, Melón, Zumos y Otras frutas. Los datos se encuentran en el archivo vegetables.zip. Aunque se pueden cargar desde el archivo vegetables.mat de Matlab, se recomienda leerlos desde los otros archivos. El archivo vegetables.dat contiene tres columnas de datos. La primera columna es la variable illness (enfermedad), la segunda es la variable veg.fruit (verduras/frutas). Y, por último, la tercera columna es la variable region (región). Los nombres de todos los elementos que necesitamos se hallan en el archivo vegetables_labels.txt de forma ordenada según la codificación.
#cargamos los datos
vegetables <- read.table("vegetables.dat")
noms <- read.table("vegetables_labels.txt")
#eliminamos los datos en blanco
rm.blanks <- function(x) gsub(" ","",x)
noms <- apply(noms,1,rm.blanks)
colnames(vegetables) <- c("illness","veg.fruit","region")
vegetables$illness <- factor(vegetables$illnes, labels=noms[1:11])
vegetables$veg.fruit <- factor(vegetables$veg.fruit, labels=noms[12:20])
vegetables$region <- factor(vegetables$region, labels=noms[21:22])
#hacemos tablas de contingencia para cada region.
attach(vegetables)
tabla.EU <- table(illness[region=="EU"],veg.fruit[region=="EU"])
tabla.EU
##
## Salad Leafy Tomato Other-vegetables Sprouts Berries
## Norovirus 15 26 1 9 0 55
## Salmonella 8 12 1 3 11 0
## E-coli 3 0 0 1 3 0
## Campylobacter 2 1 0 0 0 0
## Shigella 1 0 0 3 0 0
## Clostridium 0 1 0 6 0 0
## Staphylococcus 0 0 0 3 1 0
## Yersinia 0 0 0 3 0 0
## Bacillus 2 0 0 3 0 1
## Other-viruses 2 0 0 0 2 0
## Other-microorganisms 0 0 0 9 0 0
##
## Melon Juices Other-fruits
## Norovirus 0 0 2
## Salmonella 1 1 3
## E-coli 0 0 0
## Campylobacter 0 0 0
## Shigella 0 0 1
## Clostridium 0 0 0
## Staphylococcus 0 0 0
## Yersinia 0 0 0
## Bacillus 0 0 0
## Other-viruses 0 0 1
## Other-microorganisms 0 0 0
Tabla de contingencia para USA
tabla.USA <- table(illness[region=="USA"],veg.fruit[region=="USA"])
tabla.USA
##
## Salad Leafy Tomato Other-vegetables Sprouts Berries
## Norovirus 97 62 5 9 0 5
## Salmonella 8 8 17 3 14 2
## E-coli 10 22 0 0 4 2
## Campylobacter 4 2 1 0 0 0
## Shigella 1 2 0 0 0 0
## Clostridium 0 0 0 0 0 0
## Staphylococcus 2 0 0 0 0 0
## Yersinia 0 0 0 0 0 0
## Bacillus 1 0 0 0 0 0
## Other-viruses 0 1 1 1 0 0
## Other-microorganisms 0 0 0 0 2 0
##
## Melon Juices Other-fruits
## Norovirus 9 3 33
## Salmonella 14 0 5
## E-coli 0 6 2
## Campylobacter 1 0 1
## Shigella 0 0 0
## Clostridium 0 0 0
## Staphylococcus 0 0 0
## Yersinia 0 0 0
## Bacillus 0 0 1
## Other-viruses 0 0 2
## Other-microorganisms 1 0 0
library(ca)
ca.EU <- ca(tabla.EU)
round(ca.EU$sv^2, 3) # valores propios
## [1] 0.501 0.425 0.069 0.042 0.035 0.003 0.000 0.000
plot(ca.EU)
Tambien realizamos un plot asimetrico
plot(ca.EU, map = "rowprincipal")
Del mismo modo se puede hacer el otro gráfico asimétrico:
plot(ca.EU, map = "colprincipal")
Realizamos los mismo graficos para USA, pero antes eliminamos datos faltantes
ca.USA <- ca(tabla.USA[-c(6,8),])
round(ca.USA$sv^2, 3) # valores propios
## [1] 0.384 0.139 0.041 0.022 0.008 0.003 0.001 0.000
plot(ca.USA)
Tambien hacemos graficos asimetricos
plot(ca.USA, map = "rowprincipal")
plot(ca.USA, map = "colprincipal")
Para USA los dos primeros componentes explican el 87.2% de los datos, mientras que para la UE explican el 86.1%.
El Norovirus se muestra como responsable de la mayoría de los brotes relacionados con productos agrícolas, seguido de la Salmonela. El Norovirus está relacionado principalmente con el consumo de ensalada en los Estados Unidos y de bayas en la Unión Europea. La Salmonella fue la principal causa de brotes de productos agrícolas en varios estados de Estados Unidos y fue el patógeno implicado en la mayoría de los brotes asociados al consumo de coles, tomates y melones. Como se refleja en los gráficos, el patrón de los brotes de productos frescos difería en Estados Unidos y en la Unión Europea según el tipo de microorganismo y el vehículo alimentario implicado.
Corradino[2] utilizó el MDS para estudiar la estructura de proximidades en una colonia de monos japoneses. Las observaciones se realizaron en un grupo social de 14 monos japoneses durante un período de un año. Los catorce monos se nombran y describen en la Tabla 1. Se observaron las relaciones de proximidad cada 60 segundos. Si dos monos se encontraban a menos de 1.5 metros de distancia y se toleraban, se decía que estaban “cerca”. Se calcularon las disimilaridades para cada par de monos en función del tiempo que la pareja estuvo cerca la una de la otra. Las disimilaridades se estudiaron por separado en la temporada de cría (monk_85.dis) y fuera de ésta (monk_84.dis).
Lo mismo para el archivo monk_85.dis con las disimilaridades en la temporada de cría.
El primer paso es crear la matriz de distancia para los datos
monkeys <- read.csv("monkeys.dat", sep=";")
monkeys$age <- factor(monkeys$age)
monkeys$sex <- factor(monkeys$sex)
df <- read.table("monk_84.dis", header = FALSE, sep="", skip=2)
disim <- as.numeric(df$V1[1:91])
mat <- matrix(0, ncol=14, nrow=14)
mat[lower.tri(mat)] <- disim
mat <- mat + t(mat)
colnames(mat) <- rownames(mat) <- monkeys$alias
d.84 <- as.dist(mat)
d.84
## ALFA FRAN FELL PANC ISA GILD BETI OLGA ORSE ROSS DIVO CIST ELET
## FRAN 126
## FELL 55 211
## PANC 367 60 26
## ISA 55 39 21 12
## GILD 63 39 17 23 1
## BETI 37 74 73 32 11 3
## OLGA 27 33 24 2 28 3 3
## ORSE 96 90 68 8 14 26 24 6
## ROSS 62 30 8 51 20 1 35 8 42
## DIVO 27 18 51 100 10 17 8 5 28 0
## CIST 4 10 100 38 1 3 3 17 39 3 16
## ELET 3 19 38 16 0 35 7 16 4 4 62 6
## EVA 22 5 16 46 8 9 1 6 46 16 9 26 436
Realizamos la matriz de distancias para la otra temporada.
df <- read.table("monk_85.dis", header = FALSE, sep="", skip=2)
disim <- as.numeric(df$V1[1:91])
mat <- matrix(0, ncol=14, nrow=14)
mat[lower.tri(mat)] <- disim
mat <- mat + t(mat)
colnames(mat) <- rownames(mat) <- monkeys$alias
d.85 <- as.dist(mat)
d.85
## ALFA FRAN FELL PANC ISA GILD BETI OLGA ORSE ROSS DIVO CIST ELET
## FRAN 80
## FELL 44 161
## PANC 156 120 27
## ISA 49 18 10 21
## GILD 109 117 22 4 2
## BETI 77 113 23 40 78 11
## OLGA 83 105 9 9 7 74 18
## ORSE 53 153 105 41 25 47 69 4
## ROSS 33 44 7 11 4 2 14 16 9
## DIVO 20 29 28 75 20 44 60 14 39 0
## CIST 22 22 20 9 4 56 107 28 62 15 11
## ELET 286 7 27 7 2 41 7 14 3 15 109 34
## EVA 59 17 256 18 12 38 0 146 96 20 12 22 204
library(ade4)
is.euclid(d.84,print=F)
## Warning in is.euclid(d.84, print = F): Zero distance(s)
## [1] FALSE
#como nos faltan datos chequeamos donde estan
which(d.84==0)
## [1] 54 82
Reemplazamos estos datos por un valor pequeño cercano al "0"
d.84[which(d.84==0)] <- 0.0001
which(d.85==0)
## [1] 70 82
d.85[which(d.85==0)] <- 0.0001
Comprobamos si la distancia es euclidea
is.euclid(d.84,print=F)
## [1] FALSE
En este caso la distancia no es euclidea.
is.euclid(d.85,print=F)
## [1] FALSE
En este caso tampoco es euclidea.
Cuando la disimilaridad es euclídea o una medida directa relativa a ella (tal como el coseno), la mejor elección para el MDS es el MDS métrico. Sin embargo, cuando la disimilaridad no es euclídea o incluso no métrica, debemos decidir entre la solución MDS métrica y la no métrica. El MDS no métrico frecuentemente demuestra ser más realista y mucho mejor en la práctica.
set.seed(123)
library(MASS)
iso.mds.84 <- isoMDS(d.84)
## initial value 62.369024
## final value 62.351900
## converged
iso.mds.84$stress
## [1] 62.3519
iso.mds.85 <- isoMDS(d.85)
## initial value 47.386485
## final value 47.345903
## converged
iso.mds.85$stress
## [1] 47.3459
library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-7
nmds.84 <- metaMDS(d.84, engine="monoMDS", autotransform=FALSE, trace = FALSE)
nmds.84$stress*100
## [1] 20.91301
nmds.85 <- metaMDS(d.85, engine="monoMDS", autotransform=FALSE, trace = FALSE)
nmds.85$stress*100
## [1] 22.84126
Aunque con valores elevados de stress, es evidente que la mejor solución se obtiene con la función metaMDS().
plot(nmds.84$points[,1], nmds.84$points[,2], type = "n",
xlab="NMDS1", ylab="NMDS2")
title("Temporada sin cría")
text(nmds.84$points[,1], nmds.84$points[,2],
col=ifelse(monkeys$sex=="male","blue","red"),
labels = rownames(nmds.84$points), cex=0.5)
plot(nmds.84$points[,1], nmds.84$points[,2], type = "n",
xlab="NMDS1", ylab="NMDS2")
title("Temporada de cría")
text(nmds.85$points[,1], nmds.85$points[,2],
col=ifelse(monkeys$sex=="male","blue","red"),
labels = rownames(nmds.85$points), cex=0.5)
Observamos una mayor proximidad en la representación de la temporada de cría. Los machos aparecen agrupados en las dos configuraciones, aunque con alguna excepción. Las hembras también forman un grupo, es decir, se relacionan más entre ellas en las dos temporadas.
En el trabajo de Corradino se deducía que los machos estaban más cerca en la temporada de cría que fuera de ella. En nuestros gráficos este efecto no aparece.
¿Cual es el mono con menor residuo entre las dos configuraciones?
proc <- procrustes(nmds.84,nmds.85)
plot(proc)
El menor residuo de la representación es para Gilda.
sort(residuals(proc))[1:3]
## GILD EVA DIVO
## 31.67155 41.72386 61.20314
En primer lugar procedemos a calcular el MDS no métrico con 3 coordenadas.
nmds.85.3 <- metaMDS(d.85, engine="monoMDS", k = 3, autotransform=FALSE, trace = FALSE)
nmds.85.3$stress*100
## [1] 14.81673
El stress es más bajo que con dos dimensiones.
library(plotly)
## Loading required package: ggplot2
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:MASS':
##
## select
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
datos <- data.frame(nmds.85.3$points)
fig3D <- plot_ly(data=datos, x = ~MDS1, y = ~MDS2, z = ~MDS3,
color = ~monkeys$sex, colors = c("red","blue"))
fig3D
## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
Con la rotación adecuada se puede ver la separación total entre machos y hembras.
Normalmente, el análisis de componentes principales se lleva a cabo calculando los valores propios y los vectores propios de la matriz de correlaciones. Con n casos y p variables, si escribimos X para la matriz n×p estandarizada para que las columnas tengan media cero y desviación estándar unitaria, encontramos los valores propios y los vectores propios de la matriz p × p X0X (que es n o n − 1 veces la matriz de correlaciones dependiendo de qué denominador se haya utilizado al calcular las desviaciones estándar). Con los datos de los gorriones del ejercicio 3 de los ejercicios propuestos para PCA, comprobar las siguientes afirmaciones:
load("gorriones.RData")
X <- gorriones[,1:5]
X <- scale(X)
eig <- eigen(crossprod(X), symmetric = TRUE)
pca <- princomp(X)
pca$loadings
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## x1 0.452 0.690 0.420 0.374
## x2 0.462 -0.300 0.341 -0.548 -0.530
## x3 0.451 -0.325 -0.454 0.606 -0.343
## x4 0.471 -0.185 -0.411 -0.388 0.652
## x5 0.398 0.876 -0.178 -0.192
##
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## SS loadings 1.0 1.0 1.0 1.0 1.0
## Proportion Var 0.2 0.2 0.2 0.2 0.2
## Cumulative Var 0.2 0.4 0.6 0.8 1.0
eig$vectors
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.4517989 0.05072137 0.6904702 0.42041399 -0.3739091
## [2,] 0.4616809 -0.29956355 0.3405484 -0.54786307 0.5300805
## [3,] 0.4505416 -0.32457242 -0.4544927 0.60629605 0.3427923
## [4,] 0.4707389 -0.18468403 -0.4109350 -0.38827811 -0.6516665
## [5,] 0.3976754 0.87648935 -0.1784558 -0.06887199 0.1924341
scores <- X %*% eig$vectors
head(scores)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.06428901 -0.60083713 -0.1712334 -0.515825561 0.5487904
## [2,] -2.18031283 -0.44230082 0.4000696 -0.645459959 0.2310766
## [3,] -1.14556567 0.01925412 -0.6761269 -0.716298164 0.2088714
## [4,] -2.31106565 0.17199267 -0.3059621 0.149289289 0.4781034
## [5,] -0.29504203 -0.66520783 -0.4742138 -0.545862110 0.2444780
## [6,] 1.91626198 -0.59525444 0.6209330 0.006608669 -0.2855166
head(pca$scores)
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## [1,] 0.06428901 -0.60083713 -0.1712334 -0.515825561 -0.5487904
## [2,] -2.18031283 -0.44230082 0.4000696 -0.645459959 -0.2310766
## [3,] -1.14556567 0.01925412 -0.6761269 -0.716298164 -0.2088714
## [4,] -2.31106565 0.17199267 -0.3059621 0.149289289 -0.4781034
## [5,] -0.29504203 -0.66520783 -0.4742138 -0.545862110 -0.2444780
## [6,] 1.91626198 -0.59525444 0.6209330 0.006608669 0.2855166
D <- crossprod(scores)
eig$values
## [1] 173.566961 25.512196 18.548378 14.475145 7.897321
diag(D)
## [1] 173.566961 25.512196 18.548378 14.475145 7.897321
diag(D)/49
## [1] 3.5421829 0.5206571 0.3785383 0.2954111 0.1611698
pca$sdev^2
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## 3.5421829 0.5206571 0.3785383 0.2954111 0.1611698
load("SNP.RData")
controls <- rownames(subject.support)[subject.support$cc==0]
pop <- subject.support[controls,"stratum"]
pop.all <- subject.support[,"stratum"]
use <- seq(1, ncol(snps.10), 10)
## Loading required package: snpStats
## Loading required package: survival
## Loading required package: Matrix
ctl10 <- snps.10[controls, use]
all10 <- snps.10[,use]
X <- as(ctl10, "numeric")
X.all <- as(all10,"numeric")
dim(X)
## [1] 500 2851
La primera dificultad es substituir los valores missing por la media.
for(i in 1:ncol(X)){
X[is.na(X[,i]), i] <- mean(X[,i], na.rm = TRUE)
}
Estandarizamos los datos
X <- scale(X)
Calculamos la matriz XX'
XXt <- tcrossprod(X)
head(XXt[1:10,1:10])
## jpt.869 jpt.862 jpt.948 ceu.564 ceu.904 ceu.665 jpt.663
## jpt.869 2601.7764 138.9807 192.9046 -283.5339 -414.6340 -235.4016 247.0385
## jpt.862 138.9807 2775.1951 164.8674 -218.8962 -306.4598 -280.7771 202.4891
## jpt.948 192.9046 164.8674 2531.5733 -361.6866 -311.3615 -278.4234 309.6210
## ceu.564 -283.5339 -218.8962 -361.6866 2908.3334 439.0646 329.4606 -162.9436
## ceu.904 -414.6340 -306.4598 -311.3615 439.0646 2935.7316 313.8808 -296.6975
## ceu.665 -235.4016 -280.7771 -278.4234 329.4606 313.8808 2836.9213 -326.4090
## ceu.977 jpt.637 ceu.897
## jpt.869 -311.0834 329.9604 -95.97765
## jpt.862 -345.0864 304.8939 -306.84784
## jpt.948 -284.2715 316.8407 -261.73112
## ceu.564 260.8681 -310.9871 160.09227
## ceu.904 229.0223 -403.5789 308.65566
## ceu.665 355.8414 -296.9524 329.36431
Hacer el mismo gráfico con la segunda y tercera componentes. ¿Qué componente separa mejor las poblaciones?
eig <- eigen(XXt, symmetric = TRUE)
eig$values[1:5]
## [1] 146789.918 8589.204 8380.229 8236.379 8111.603
U <- eig$vectors
oldpar <- par(mfrow=c(1,3))
boxplot(U[,1] ~ pop, ylab="PC1")
boxplot(U[,2] ~ pop, ylab="PC2")
boxplot(U[,3] ~ pop, ylab="PC3")
par(oldpar)
eig$values[eig$values<0] <- 0
v <- 1/sqrt(eig$values)
m <- crossprod(X, U)
# B <- m %*% diag(v) # poco eficiente
# B <- t(t(m) * v) # más eficiente
B <- m * rep(v, rep(nrow(m), ncol(m))) # mejor
B[1:10,1:10]
## [,1] [,2] [,3] [,4] [,5]
## rs7909677 0.0003252226 -9.331756e-03 0.015445805 0.003947906 -0.011546724
## rs4880781 -0.0055876444 7.507270e-03 -0.000543355 -0.008346595 -0.003828963
## rs2448365 -0.0076646829 1.733169e-02 -0.008138940 0.007698970 0.007641585
## rs7919436 0.0209014311 1.102595e-02 -0.025856734 0.013353554 0.024974745
## rs12357593 0.0301087215 1.708045e-03 -0.022693905 0.014494197 0.022537829
## rs10904173 -0.0032615624 1.569638e-02 0.017220464 0.002238865 -0.014006605
## rs11252693 0.0107882727 -4.632472e-05 -0.012653201 0.011259890 0.006176045
## rs11253096 0.0259380458 -1.762775e-02 0.006504495 0.010741717 -0.010751181
## rs17159711 0.0150469759 1.284841e-02 0.010802802 0.007726009 0.030637888
## rs4881518 0.0168951870 -1.775938e-02 -0.008888053 0.018592907 0.004927426
## [,6] [,7] [,8] [,9] [,10]
## rs7909677 0.009234650 0.009592387 -0.0006004881 0.009523450 -0.014859529
## rs4880781 -0.025199111 0.001262191 0.0136031063 -0.034072213 -0.011145068
## rs2448365 -0.009577234 0.025332822 0.0176713475 0.005336778 -0.009796594
## rs7919436 0.004708050 -0.017783670 0.0129060916 0.035684897 -0.009938182
## rs12357593 0.017645674 -0.035232563 0.0028693908 0.036582396 -0.002629847
## rs10904173 0.008486269 0.013840947 0.0097519493 -0.023745230 -0.017696460
## rs11252693 0.019020146 -0.019180868 -0.0004079738 0.033007216 0.014243545
## rs11253096 0.004788495 -0.014961988 0.0022505070 -0.024822305 -0.008882354
## rs17159711 -0.011462165 0.022790028 -0.0054857879 0.008332284 -0.009762654
## rs4881518 0.022255034 -0.006821502 0.0213503830 0.025230784 0.001849588
Se controlan los valores faltantes y se estandarizan los datos
for(i in 1:ncol(X.all)){
X.all[is.na(X.all[,i]), i] <- mean(X.all[,i], na.rm = TRUE)
}
X.all <- scale(X.all)
Calculamos las puntuaciones y hacemos el grafico
scores.all.12 <- X.all %*% B[,1:2]
plot(scores.all.12, xlab="PC1", ylab="PC2",
col=ifelse(pop.all=="CEU","red","blue"))
Las dos poblaciones quedan también separadas por la primera componente.